﻿using System;
using System.IO;

namespace Automata.examen2
{
    class EigenVectores
    {
        private double[,] matriz;
        private double[,] balancedMatrix;
        private double[,] hessenbergMatrix;
        private double[] aprox;
        private double[] vector;

        private int N;
        private double ultimoError;
        private int totalIteraciones = 0;

        private const int MAX_ITERACIONES = 500;
        private const double LIMITE = 0.001;

        private int totalMultiplicaciones = 0;

        private Complex[] eigenValuesComplex;
        private double[] eigenValues;

        /// <summary>
        /// Carga la matriz del archivo en la viable matriz
        /// </summary>
        /// <param name="archivo">El stream del archivo donde se encuentra las matriz</param>
        private void cargaMatriz(Stream archivo)
        {
            string contenido = EntradaSalida.leeArchivo(archivo);

            string[] lineas = contenido.Split(new char[] { ';' });

            N = lineas.Length;

            matriz = new double[N, N];
            aprox = new double[N];
            vector = new double[N];

            for (int i = 0; i < lineas.Length; i++)
            {
                string[] celdas = lineas[i].Split(new char[] { ' ' });

                for (int j = 0; j < lineas.Length; j++)
                    matriz[i, j] = Double.Parse(celdas[j]);

                aprox[i] = 1;
            }
        }

        /// <summary>
        /// Calcula el eigenvalor y eigenvector dominante de la matriz si existe
        /// </summary>
        /// <param name="archivo">El archivo con la matriz</param>
        public string genera(Stream archivo)
        {
            try
            {
                cargaMatriz(archivo);
            }
            catch (Exception)
            {
                return "El archivo no tiene el formato correcto";
            }

            return aproxima();
        }

        /// <summary>
        /// Aproxima los eigenvalores e eigenvecores
        /// </summary>
        /// <returns>Un string amigable con los eigenvectores e igenvalores</returns>
        private string aproxima()
        {
            string vec = eigenvector();
            string val = eigenvalor();

            return "Con el método de <i>Power Iteration</i> se obtuvo el eigenvector dominante <br /> <br /> " + vec +
                   "<br /><br /> y el eigenvalor dominante: <br /> <br /> " + val +
                   "<br /><br /> Usando " + totalMultiplicaciones + " multiplicaciones, con un error de " + ultimoError +
                   " y " + totalIteraciones + " iteraciones.";
        }

        /// <summary>
        /// Multiplica las matrices por fuerza bruta y devuelve la matriz
        /// aproximanda de la iteracion
        /// </summary>
        /// <returns>La matriz multiplicada</returns>
        private double[] multiplica()
        {
            totalMultiplicaciones = 0;

            double[] temp = new double[N];

            for (int i = 0; i < N; i++)
                for (int k = 0; k < N; k++)
                {
                    if (matriz[i, k] == 0 || aprox[k] == 0)
                        continue;
                    temp[i] += matriz[i, k] * aprox[k];
                    totalMultiplicaciones++;
                }

            return temp;
        }

        /// <summary>
        /// Estandariza los valores de temp y compara la diferencia con aprox
        /// si es menor que el límite, se dice que ya acabo
        /// </summary>
        /// <returns>Si ya acabamos o no</returns>
        private bool yaAcabamos(double[] temp)
        {
            double menor = temp[0];
            double error = 0;
            double[] act = new double[N];

            foreach (long i in temp)
            {
                if (Math.Abs(i) < Math.Abs(menor))
                    menor = i;
            }

            for (int i = 0; i < N; i++)
            {
                act[i] = temp[i] / menor;
                error += Math.Abs(act[i] - vector[i]);
            }

            vector = act;
            ultimoError = error;

            return (error < LIMITE);
        }

        /// <summary>
        /// Aproxima el eigenvector
        /// </summary>
        /// <returns>Un string con el eigenvector</returns>
        private string eigenvector()
        {
            int iteraciones = 0;
            double[] temp = null;

            while (iteraciones < MAX_ITERACIONES)
            {
                temp = multiplica();
                if (yaAcabamos(temp))
                    break;
                aprox = temp;
                iteraciones++;
                totalIteraciones++;
            }

            if (iteraciones == MAX_ITERACIONES)
                return "La matriz introducida no tiene un eigenvalor dominante por lo que este método no puede encontrar ningún eigenvalor.";

            string s = "[";

            for (int i = 0; i < (N - 1); i++)
                s += vector[i] + ", ";

            s += vector[N - 1] + "]";

            return s;
        }

        /// <summary>
        /// Calcula el producto punto de m1 y m2
        /// </summary>
        /// <param name="m1">La matriz m1</param>
        /// <param name="m2">La matriz m2</param>
        /// <returns>El producto punto</returns>
        private double punto(double[] m1, double[] m2)
        {
            double tot = 0;

            for (int i = 0; i < N; i++)
            {
                tot += m1[i] * m2[i];
                totalMultiplicaciones++;
            }

            return tot;
        }

        /// <summary>
        /// Calcula el eigenvalor del eigenvector calculado
        /// </summary>
        /// <returns>El eigenvalor</returns>
        private string eigenvalor()
        {
            aprox = vector;

            double valor = punto(multiplica(), vector) / punto(vector, vector);

            return valor.ToString();
        }

        public void balanceMatrix()
        {
            int RADIX = 2;
            int i, j, last = 0;
            double s, r, g, f, c, sqrdx;

            double[,] matrixCopy=this.matriz;

            sqrdx = RADIX * RADIX;
            while (last == 0)
            {
                last = 1;
                for (i = 0; i < N; i++)
                {
                    r = c = 0.0;
                    for (j = 0; j < N; j++)
                    {
                        if (j != i)
                        {
                            c+=Math.Abs(matrixCopy[j,i]);
                            r += Math.Abs(matrixCopy[i, j]);
                        }
                    }
                    if(c!=0.0 && r!= 0.0)
                    {
                        g = r / RADIX;
                        f = 1.0;
                        s = c + r;
                        while (c < g)
                        {
                            f *= RADIX;
                            c *= sqrdx;
                        }
                        g = r * RADIX;
                        while (c > g)
                        {
                            f /= RADIX;
                            c /= RADIX;
                        }
                        if ((c + r) / f < 0.95 * s)
                        {
                            last = 0;
                            g = 1.0 / f;
                            for (j = 0; j < N; j++) matrixCopy[i, j] *= g;
                            for (j = 0; j < N; j++) matrixCopy[j, i] *= f;
                        }
                    }
                }
            }
            this.balancedMatrix = matrixCopy;
        }

        /// <summary>
        /// Calcula la matriz de hessenberg the la matriz actual y lo asigna a la variable hessenbergMatrix
        /// </summary>
        public void createHessenberg()
        {
            double[,] matrixCopy = this.balancedMatrix;
            hessenbergMatrix = new double[N, N];

            // numerical recipes
            int i, j, m;
            double y, x;

            for (m = 1; m < N; m++)
            {
                x = 0.0;
                i = m;
                //Find the pivot
                for (j = m; j < N; j++)
                {
                    if (Math.Abs(matrixCopy[j, m - 1]) > Math.Abs(x))
                    {
                        x = matrixCopy[j, m - 1];
                        i = j;
                    }
                }
                //Interchange rows and columns
                if (i != m)
                {
                    for (j = m - 1; j < N; j++)
                    {
                        y = matrixCopy[i, j];
                        matrixCopy[i, j] = matrixCopy[m, j];
                        matrixCopy[m, j] = y;
                    }
                    for (j = 0; j < N; j++)
                    {
                        y = matrixCopy[j, i];
                        matrixCopy[j, i] = matrixCopy[j, m];
                        matrixCopy[j, m] = y;
                    }
                }
                //Carry out the elimination
                if (x != 0.0)
                {
                    for (i = m + 1; i < N; i++)
                    {
                        if ((y = matrixCopy[i, m - 1]) != 0.0)
                        {
                            y /= x;
                            matrixCopy[i, m - 1] = y;
                            for (j = m; j < N; j++)
                            {
                                matrixCopy[i, j] -= y * matrixCopy[m, j];
                            }
                            for (j = 0; j < N; j++)
                            {
                                matrixCopy[j, m] += y * matrixCopy[j, i];
                            }
                        }
                    }
                }
            }
            this.hessenbergMatrix = matrixCopy;
        }//end hessenberg
        

        private double Sign(double a, double b)
        {
            return ((b) >= 0.0 ? Math.Abs(a) : -Math.Abs(a));
        }

        public void hqr()
        {
            int nn, m, l, k, j, its, i, mmin;
            double z,y,x,w,v,u,t,s,r,q,p,anorm;
            double [,] hessenbergCopy=hessenbergMatrix;
            double[] wri = new double[N];
            Complex[]wric=new Complex[N];

            anorm = 0.0;
            p = 0.0;
            q = 0.0;
            r = 0.0;
            //Compute matrix norm for possible us in locating
            //single small subdiagonal elements
            for (i = 0; i < N; i++)
            {
                for (j = Math.Max(i - 1, 0); j < N; j++)
                {
                    anorm += Math.Abs(hessenbergCopy[i, j]);
                }
            }
            nn = N - 1;
            t = 0.0;
            //Begin search for next eigenvalue
            while (nn >= 0)
            {
                its = 0;
                do{
                    for (l = nn; l > 0; l--)
                    {
                        s = Math.Abs(hessenbergCopy[l-1,l-1]+Math.Abs(hessenbergCopy[l,l]));
                        if (s == 0.0) s = anorm;
                        if(Math.Abs(hessenbergCopy[l,l-1]) + s ==s)
                        {
                            hessenbergCopy[l, l - 1] = 0.0;
                            break;
                        }
                    }
                    x = hessenbergCopy[nn,nn];
                    //One root found
                    if (l == nn)
                    {
                        wri[nn--] = x + t;
                    }
                    else
                    {
                        y = hessenbergCopy[nn - 1, nn - 1];
                        w = hessenbergCopy[nn, nn - 1] * hessenbergCopy[nn - 1, nn];
                        //Two roots found
                        if (l == nn - 1)
                        {
                            p = 0.5 * (y - x);
                            q = p * p + w;
                            z = Math.Sqrt(Math.Abs(q));
                            x += t;
                            //a real pair
                            if (q >= 0.0)
                            {
                                z = p + Sign(z, p);
                                wri[nn - 1] = wri[nn] = x + z;
                                if (z != 0.0) wri[nn] = x - w / z;
                            }
                            //a complex pair
                            else
                            {
                               wric[nn] = new Complex(x + p, -z);
                               wric[nn - 1] = Complex.conjg(wric[nn]);
                            }
                            nn -= 2;
                        }
                        //No roots found. Continue Iteration
                        else
                        {
                            if (its == 30) throw new Exception("Too many iterations in hqr");
                            if (its == 10 || its == 20)
                            {
                                t += x;
                                for (i = 0; i < nn + 1; i++) hessenbergCopy[i, i] -= x;
                                s=Math.Abs(hessenbergCopy[nn,nn-1]+Math.Abs(hessenbergCopy[nn-1,nn-2]));
                                y = x = 0.75 * s;
                                w = -0.4375 * s * s;
                            }
                            ++its;
                            for(m=nn-2;m>=1;m--){
                                z=hessenbergCopy[m,m];
                                r=x-z;
                                s=y-z;
                                p=(r*s-w)/hessenbergCopy[m+1,m]+hessenbergCopy[m,m+1];
                                q=hessenbergCopy[m+1,m+1]-z-r-s;
                                r=hessenbergCopy[m+2,m+1];
                                s = Math.Abs(p) + Math.Abs(q) + Math.Abs(r);
                                p /= s;
                                q /= s;
                                r /= s;
                                if (m == 1) break;
                                u = Math.Abs(hessenbergCopy[m, m - 1]) * (Math.Abs(q) + Math.Abs(r));
                                v = Math.Abs(p)*(Math.Abs(hessenbergCopy[m-1,m-1]+Math.Abs(z)+Math.Abs(hessenbergCopy[m+1,m+1])));
                                if(u+v ==v) break;
                            }
                            for(i=m;i<nn-1;i++)
                            {
                                hessenbergCopy[i+2,i]=0.0;
                                if(i !=m) hessenbergCopy[i+2,i+1]=0.0;
                            }
                            //Double QR step on rows 1 to nn and columns m to nn
                            for(k=m;k<nn;k++)
                            {
                                if(k !=m)
                                {
                                    p=hessenbergCopy[k,k-1];
                                    q=hessenbergCopy[k+1,k+1];
                                    r=0.0;
                                    if(k+1 != nn) r=hessenbergCopy[k+2,k-1];
                                    if((x=Math.Abs(p)+Math.Abs(q)+Math.Abs(r))!=0.0)
                                    {
                                        p /= x;
                                        q /= x;
                                        r /= x;
                                    }
                                }
                                if((s=Sign(Math.Sqrt(p*p+q*q+r*r),p))!=0.0)
                                {
                                    if (k == m)
                                    {
                                        if (l != m) hessenbergCopy[k, k - 1] = -hessenbergCopy[k, k - 1];
                                    }
                                    else
                                    {
                                        hessenbergCopy[k, k - 1] = -s * x;
                                    }
                                    p += s;
                                    x = p / s;
                                    y = q / s;
                                    z = r / s;
                                    q /= p;
                                    r /= p;
                                    
                                    for(j=k;j<nn+1;j++)
                                    {
                                        p=hessenbergCopy[k,j]+q*hessenbergCopy[k+1,j];
                                        if (k+1 != nn) {
                                            p += r*hessenbergCopy[k+2,j];
                                            hessenbergCopy[k+2,j] -= p*z;
                                        }
                                        hessenbergCopy[k+1,j] -= p*y;
                                        hessenbergCopy[k,j] -= p*x;
                                    }
                                    mmin = nn < k+3 ? nn : k+3;
                                    for (i=l;i<mmin+1;i++) 
                                    {
                                        p=x*hessenbergCopy[i,k]+y*hessenbergCopy[i,k+1];
                                        if (k+1 != (nn)) 
                                        {
                                            p += z*hessenbergCopy[i,k+2];
                                            hessenbergCopy[i,k+2] -= p*r;
                                        }
                                        hessenbergCopy[i,k+1] -= p*q;
                                        hessenbergCopy[i,k] -= p;
                                    }
                                }
                            }
                        }
                    }
                }while(l+1<nn);
            }
            this.eigenValues=wri;
            this.eigenValuesComplex=wric;
        }//end hqr
    }
}
